Chapter 5 Community composition

load("data/data.Rdata")

5.1 Taxonomy overview

5.1.1 Stacked barplot

# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  filter(count > 0) #filter 0 counts

# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)

ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(. ~ type,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")+
  scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
  facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label

5.1.2 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_96tlbr7yta7x19qx594q
phylum mean sd
p__Bacteroidota 0.383073524 0.200556253
p__Bacillota_A 0.255324749 0.163209517
p__Bacillota 0.118782585 0.146832611
p__Pseudomonadota 0.094122159 0.156932597
p__Campylobacterota 0.053366123 0.093185234
p__Verrucomicrobiota 0.027449290 0.066220962
p__Desulfobacterota 0.023720834 0.036545754
p__Chlamydiota 0.010635637 0.059863092
p__Fusobacteriota 0.010544526 0.028384721
p__Cyanobacteriota 0.009261265 0.016525997
p__Bacillota_C 0.004741218 0.006655526
p__Spirochaetota 0.004033548 0.012340900
p__Bacillota_B 0.002480140 0.004938693
p__Actinomycetota 0.001243097 0.006365104
p__Elusimicrobiota 0.001221306 0.006520508
phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    select(phylum) %>%
    pull()

phylum_summary %>%
    filter(phylum %in% phylum_arrange) %>%
    mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
    ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
        scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        geom_jitter(alpha=0.5) + 
        theme_minimal() + 
        theme(legend.position="none") +
        labs(y="Phylum",x="Relative abundance")

5.2 Taxonomy boxplot

5.2.1 Family

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_summary %>%
    group_by(family) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_5sgigi9whnyw943pevpo
family mean sd
f__Bacteroidaceae 2.228331e-01 0.1385585120
f__Lachnospiraceae 1.414562e-01 0.1083170968
f__Tannerellaceae 1.034197e-01 0.0797797418
f__Helicobacteraceae 5.290294e-02 0.0927590950
f__Mycoplasmoidaceae 3.716482e-02 0.0758142978
f__Erysipelotrichaceae 3.523239e-02 0.0451748821
f__UBA3700 3.438944e-02 0.0558522623
f__Marinifilaceae 2.790560e-02 0.0271961340
f__DTU072 2.716390e-02 0.0530196982
f__Rikenellaceae 2.712953e-02 0.0465396092
f__Enterobacteriaceae 2.707699e-02 0.0920895535
f__Coprobacillaceae 2.566962e-02 0.0894627842
f__ 2.527324e-02 0.0777603339
f__Desulfovibrionaceae 2.372083e-02 0.0365457537
f__Ruminococcaceae 1.871009e-02 0.0430525690
f__UBA3830 1.588214e-02 0.0455209161
f__Rhizobiaceae 1.532665e-02 0.0768399323
f__LL51 1.510265e-02 0.0608228261
f__Akkermansiaceae 1.234664e-02 0.0317643089
f__Chlamydiaceae 1.063564e-02 0.0598630919
f__Fusobacteriaceae 1.054453e-02 0.0283847213
f__CAG-239 9.195743e-03 0.0151277476
f__Enterococcaceae 8.470645e-03 0.0467930117
f__Gastranaerophilaceae 7.774775e-03 0.0144711894
f__Oscillospiraceae 6.682761e-03 0.0078169630
f__UBA1997 6.416374e-03 0.0310554945
f__Streptococcaceae 6.404337e-03 0.0343162283
f__UBA1242 4.701177e-03 0.0154146977
f__Brevinemataceae 4.033548e-03 0.0123409004
f__Acutalibacteraceae 3.394636e-03 0.0109856431
f__UBA660 3.215538e-03 0.0139953943
f__Clostridiaceae 2.859123e-03 0.0171837280
f__RUG11792 2.834502e-03 0.0251869021
f__Peptococcaceae 2.480140e-03 0.0049386928
f__MGBC116941 2.173027e-03 0.0094092981
f__Acidaminococcaceae 1.955375e-03 0.0050449819
f__CAG-508 1.818740e-03 0.0064351937
f__Anaerovoracaceae 1.578874e-03 0.0036560185
f__Moraxellaceae 1.494257e-03 0.0097731073
f__RUG14156 1.486490e-03 0.0045970591
f__Staphylococcaceae 1.369815e-03 0.0051004179
f__Elusimicrobiaceae 1.221306e-03 0.0065205081
f__CAG-288 9.547359e-04 0.0060374248
f__Anaerotignaceae 9.043255e-04 0.0040584492
f__CALVMC01 7.561589e-04 0.0043736226
f__Eggerthellaceae 6.446025e-04 0.0021324245
f__Massilibacillaceae 6.272186e-04 0.0016392219
f__Mycobacteriaceae 5.984944e-04 0.0060578862
f__UBA1820 4.764220e-04 0.0013091229
f__Arcobacteraceae 4.631865e-04 0.0050473402
f__CAG-274 4.546649e-04 0.0022091547
f__Burkholderiaceae_C 3.721451e-04 0.0048235514
f__Muribaculaceae 3.639429e-04 0.0009801614
f__UBA932 3.439904e-04 0.0011614775
f__Hepatoplasmataceae 3.006899e-04 0.0038973865
f__Rhodobacteraceae 2.976706e-04 0.0038582514
f__Weeksellaceae 2.788124e-04 0.0031569776
f__Eubacteriaceae 1.656625e-04 0.0006747973
f__Sphingobacteriaceae 1.514738e-04 0.0012496721
f__Devosiaceae 1.498864e-04 0.0015139002
f__Pumilibacteraceae 1.285021e-04 0.0007668974
f__WRAU01 9.660522e-05 0.0012521468
f__Peptostreptococcaceae 2.300953e-05 0.0002982376
family_arrange <- family_summary %>%
    group_by(family) %>%
    summarise(mean=sum(relabun)) %>%
    arrange(-mean) %>%
    select(family) %>%
    pull()

family_summary %>%
    left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    filter(family %in% family_arrange[1:20]) %>%
    mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

5.2.2 Genus

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_summary %>%
    group_by(genus) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_z0fl9hamrbeet34h5bar
genus mean sd
g__Bacteroides 1.360888e-01 0.0920042645
g__Parabacteroides 9.738932e-02 0.0800669017
g__Phocaeicola 6.976407e-02 0.0795272001
g__Mycoplasmoides 3.069452e-02 0.0755136832
g__Helicobacter_J 3.026915e-02 0.0596129112
g__Odoribacter 2.567500e-02 0.0267638614
g__Roseburia 2.398913e-02 0.0560807126
g__NHYM01 2.263379e-02 0.0799142802
g__Alistipes 2.221496e-02 0.0283922644
g__Coprobacillus 2.025846e-02 0.0881340464
g__Agrobacterium 1.532665e-02 0.0768399323
g__Akkermansia 1.234664e-02 0.0317643089
g__Fusobacterium_A 1.045087e-02 0.0283902792
g__Kineothrix 8.853496e-03 0.0410203841
g__Proteus 8.709520e-03 0.0683836665
g__Dielma 8.629168e-03 0.0090864942
g__CAG-95 8.157620e-03 0.0205536177
g__JAAYNV01 8.042189e-03 0.0195890923
g__UBA866 7.303541e-03 0.0294323192
g__Desulfovibrio 7.059930e-03 0.0212042084
g__Enterococcus 7.042975e-03 0.0457933183
g__Ureaplasma 6.470299e-03 0.0138334733
g__Lactococcus 6.404337e-03 0.0343162283
g__Lacrimispora 6.100861e-03 0.0098371428
g__Parabacteroides_B 6.030362e-03 0.0100365790
g__CALXRO01 5.800048e-03 0.0309414387
g__Citrobacter 5.751065e-03 0.0335518785
g__NSJ-61 5.593207e-03 0.0199387071
g__Breznakia 5.537293e-03 0.0237686728
g__Clostridium_AQ 5.357893e-03 0.0121988801
g__Salmonella 5.164140e-03 0.0185065215
g__Bilophila 5.013506e-03 0.0088636099
g__Hungatella_A 4.840444e-03 0.0095762431
g__MGBC136627 4.390777e-03 0.0163455561
g__Escherichia 4.213296e-03 0.0266876301
g__UMGS1251 4.184603e-03 0.0072862642
g__Hungatella 4.140890e-03 0.0191331978
g__Clostridium_Q 4.035877e-03 0.0052190471
g__Brevinema 4.033548e-03 0.0123409004
g__Thomasclavelia 3.925809e-03 0.0109326053
g__Scatousia 3.671667e-03 0.0102995947
g__Enterocloster 3.641933e-03 0.0047397590
g__Mailhella 3.633580e-03 0.0102738868
g__Copromonas 3.620793e-03 0.0050252903
g__Ventrimonas 3.534268e-03 0.0071394041
g__Caccovivens 3.363217e-03 0.0122532776
g__Lawsonia 3.308592e-03 0.0117503555
g__Fournierella 3.236376e-03 0.0062446298
g__Limenecus 3.182108e-03 0.0065882600
g__MGBC133411 3.060527e-03 0.0074762303
g__Mucinivorans 2.917358e-03 0.0374302862
g__Sarcina 2.859123e-03 0.0171837280
g__Acetatifactor 2.754437e-03 0.0058585690
g__Eisenbergiella 2.713468e-03 0.0068503944
g__Bacteroides_G 2.698691e-03 0.0347178966
g__CAJLXD01 2.649496e-03 0.0087733544
g__Blautia 2.573938e-03 0.0061559204
g__C-19 2.289060e-03 0.0048805592
g__Butyricimonas 2.230601e-03 0.0050201300
g__Velocimicrobium 2.214326e-03 0.0066941819
g__CAZU01 2.209459e-03 0.0066119796
g__MGBC116941 2.173027e-03 0.0094092981
g__Intestinimonas 2.094290e-03 0.0039467287
g__Negativibacillus 2.081393e-03 0.0055279017
g__Rikenella 1.997207e-03 0.0037146132
g__Phascolarctobacterium 1.955375e-03 0.0050449819
g__RGIG6463 1.800497e-03 0.0039777189
g__JALFVM01 1.752732e-03 0.0038715710
g__Oscillibacter 1.519449e-03 0.0025039951
g__Pseudoflavonifractor 1.511922e-03 0.0027764599
g__Acinetobacter 1.494257e-03 0.0097731073
g__Citrobacter_A 1.401214e-03 0.0060519635
g__Staphylococcus 1.369815e-03 0.0051004179
g__RGIG4733 1.311551e-03 0.0040589196
g__UBA1436 1.221306e-03 0.0065205081
g__Lachnotalea 1.215388e-03 0.0049307664
g__Ruthenibacterium 1.209176e-03 0.0053784424
g__14-2 1.191694e-03 0.0096582871
g__Beduini 1.180938e-03 0.0025200906
g__Scatocola 1.127866e-03 0.0045101626
g__Faecisoma 1.092047e-03 0.0055756491
g__Enterococcus_A 1.090444e-03 0.0099443378
g__Faecimonas 9.939338e-04 0.0083324084
g__RGIG9287 9.696200e-04 0.0093504118
g__CAG-345 9.547359e-04 0.0060374248
g__Blautia_A 9.262837e-04 0.0029160490
g__RGIG1896 8.393274e-04 0.0062956893
g__CAG-269 8.029962e-04 0.0047313823
g__Marseille-P3106 7.985520e-04 0.0017660762
g__WRHT01 6.467835e-04 0.0027055875
g__Eggerthella 6.446025e-04 0.0021324245
g__IOR16 6.437031e-04 0.0020707005
g__Anaerotruncus 6.302739e-04 0.0016992229
g__RUG14156 6.188398e-04 0.0022208214
g__CHH4-2 6.181620e-04 0.0020051727
g__Corynebacterium 5.984944e-04 0.0060578862
g__Serratia_A 5.895501e-04 0.0076414423
g__CAG-56 4.944690e-04 0.0016386388
g__Merdimorpha 4.764220e-04 0.0013091229
g__MGBC140009 4.707187e-04 0.0024136538
g__CALURL01 4.662877e-04 0.0016783668
g__Aliarcobacter 4.631865e-04 0.0050473402
g__Scatenecus 4.547121e-04 0.0019816847
g__RGIG8482 4.425249e-04 0.0029840978
g__Enterobacter 4.097683e-04 0.0041440045
g__Klebsiella 4.078573e-04 0.0049056069
g__Caccenecus 3.964657e-04 0.0017852972
g__Alcaligenes 3.721451e-04 0.0048235514
g__Plesiomonas 3.654875e-04 0.0027184625
g__HGM05232 3.639429e-04 0.0009801614
g__JAHHSE01 3.613502e-04 0.0014942841
g__Egerieousia 3.439904e-04 0.0011614775
g__Emergencia 3.433949e-04 0.0017500504
g__Enterococcus_B 3.372270e-04 0.0022331962
g__Stoquefichus 3.044085e-04 0.0020563925
g__Hepatoplasma 3.006899e-04 0.0038973865
g__Paracoccus 2.976706e-04 0.0038582514
g__Moheibacter 2.788124e-04 0.0031569776
g__Scatomorpha 2.656736e-04 0.0010212728
g__UBA7185 2.448818e-04 0.0014600491
g__Eubacterium 1.656625e-04 0.0006747973
g__Sphingobacterium 1.514738e-04 0.0012496721
g__Devosia 1.498864e-04 0.0015139002
g__Anaerosporobacter 1.462768e-04 0.0012784872
g__Caccomorpha 1.391355e-04 0.0010571570
g__UBA2658 1.315233e-04 0.0007225792
g__Protoclostridium 1.285021e-04 0.0007668974
g__Angelakisella 1.276029e-04 0.0009248545
g__Cetobacterium_A 9.365634e-05 0.0008744341
g__Rahnella 6.509221e-05 0.0008436915
g__Peptostreptococcus 2.300953e-05 0.0002982376
genus_arrange <- genus_summary %>%
    group_by(genus) %>%
    summarise(mean=sum(relabun)) %>%
    filter(genus != "g__")%>%
    arrange(-mean) %>%
    select(genus) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    pull()

genus_summary %>%
    left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    filter(genus %in% genus_arrange[1:20]) %>%
    mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
        scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

5.3 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

5.3.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric ) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

5.3.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

5.3.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

5.3.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

5.3.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

5.3.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

5.3.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

5.4 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

5.5 Permanovas

5.5.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

5.5.1 1. Are the wild populations similar?

5.5.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

5.5.1.2 Number of samples used

[1] 28
beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

5.5.1.3 Richness


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.002848 0.0028483 0.2723    999  0.624
Residuals 26 0.272017 0.0104622                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.619
Podarcis_muralis            0.60624                 
Df SumOfSqs R2 F Pr(>F)
species 1 1.631953 0.207198 6.795072 0.001
Residual 26 6.244347 0.792802 NA NA
Total 27 7.876300 1.000000 NA NA

5.5.1.4 Neutral


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.002893 0.0028935 0.2628    999  0.598
Residuals 26 0.286298 0.0110115                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.593
Podarcis_muralis            0.61255                 
Df SumOfSqs R2 F Pr(>F)
species 1 2.018797 0.2566342 8.976049 0.001
Residual 26 5.847641 0.7433658 NA NA
Total 27 7.866438 1.0000000 NA NA

5.5.1.5 Phylogenetic


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07936 0.079362 4.7909    999  0.037 *
Residuals 26 0.43070 0.016565                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.033
Podarcis_muralis           0.037786                 
Df SumOfSqs R2 F Pr(>F)
species 1 0.3786052 0.2108638 6.947419 0.001
Residual 26 1.4168908 0.7891362 NA NA
Total 27 1.7954960 1.0000000 NA NA

5.5.1.6 Functional


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.016513 0.016513 1.4504    999  0.254
Residuals 26 0.296006 0.011385                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.277
Podarcis_muralis            0.23931                 
Df SumOfSqs R2 F Pr(>F)
species 1 0.0787916 0.1594096 4.930642 0.072
Residual 26 0.4154800 0.8405904 NA NA
Total 27 0.4942716 1.0000000 NA NA
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

5.5.2 2. Do the antibiotics work?

5.5.2.1 Acclimation vs antibiotics

treat <- meta  %>% #antibiotics samples giving error when plotting
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

5.5.2.2 Number of samples used

[1] 52
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)
5.5.2.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.027648 0.0276478 6.8746    999   0.01 **
Residuals 50 0.201088 0.0040218                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.015
2_Antibiotics      0.011555              
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.952865 0.0941315 6.574358 0.001
species 1 2.299135 0.1108223 7.740081 0.001
individual 26 9.662166 0.4657330 1.251072 0.003
Residual 23 6.831983 0.3293133 NA NA
Total 51 20.746150 1.0000000 NA NA
5.5.2.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.043782 0.043782 8.3826    999  0.004 **
Residuals 50 0.261147 0.005223                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.005
2_Antibiotics     0.0056034              
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.078576 0.1048012 9.380227 0.001
species 1 3.046369 0.1535971 13.747699 0.001
individual 26 9.611970 0.4846327 1.668348 0.001
Residual 23 5.096598 0.2569690 NA NA
Total 51 19.833512 1.0000000 NA NA
5.5.2.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.60296 0.60296 37.862    999  0.001 ***
Residuals 50 0.79626 0.01593                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.001
2_Antibiotics    1.2645e-07              
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.8994052 0.2145480 20.282274 0.001
species 1 0.8219541 0.0928441 8.777010 0.001
individual 26 3.9777800 0.4493115 1.633678 0.003
Residual 23 2.1539163 0.2432964 NA NA
Total 51 8.8530557 1.0000000 NA NA
5.5.2.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.18583 0.185833 5.3501    999   0.03 *
Residuals 50 1.73673 0.034735                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.032
2_Antibiotics      0.024874              
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.935980 0.3884761 38.1894171 0.001
species 1 0.003156 0.0006333 0.0622557 0.809
individual 26 1.878423 0.3769266 1.4251551 0.172
Residual 23 1.165966 0.2339640 NA NA
Total 51 4.983525 1.0000000 NA NA
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

5.5.3 3. Do the FMT work?

5.5.3.1 Comparison between FMT2 vs Post-FMT2

transplant3 <- meta  %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))

transplant3_nmds <- sample_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

5.5.3.2 Number of samples used

[1] 53
beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts, q=1, dist=dist)
5.5.3.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.9809525 0.0778442 5.775712 0.001
time_point 1 0.7092107 0.0562799 4.175734 0.001
type 1 1.4479860 0.1149060 8.525540 0.001
individual 25 6.9157236 0.5488022 1.628753 0.001
Residual 15 2.5476146 0.2021678 NA NA
Total 43 12.6014875 1.0000000 NA NA
5.5.3.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 1.1101711 0.0895457 7.879642 0.001
time_point 1 0.7363935 0.0593971 5.226687 0.001
type 1 1.9027421 0.1534740 13.505059 0.001
individual 25 6.5351386 0.5271203 1.855374 0.001
Residual 15 2.1133659 0.1704628 NA NA
Total 43 12.3978112 1.0000000 NA NA
5.5.3.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.1341431 0.1005265 6.921553 0.001
time_point 1 0.1159136 0.0868654 5.980943 0.001
type 1 0.1423573 0.1066822 7.345390 0.001
individual 25 0.6512838 0.4880704 1.344205 0.067
Residual 15 0.2907075 0.2178554 NA NA
Total 43 1.3344053 1.0000000 NA NA
5.5.3.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 -0.0031096 -0.0030953 -0.1177247 0.824
time_point 1 -0.0045864 -0.0045653 -0.1736359 0.833
type 1 0.0775554 0.0771987 2.9361453 0.141
individual 25 0.5385511 0.5360739 0.8155530 0.621
Residual 15 0.3962105 0.3943880 NA NA
Total 43 1.0046211 1.0000000 NA NA
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

5.5.4 4. Are there differences between the control and the treatment group?

5.5.4.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

5.5.4.2 Number of samples used

[1] 27
beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)
5.5.4.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.015090 0.0075452 2.1059    999  0.136
Residuals 24 0.085988 0.0035828                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
             Control Hot_control Treatment
Control                 0.004000     0.638
Hot_control 0.009502                 0.220
Treatment   0.624826    0.248542          
Df SumOfSqs R2 F Pr(>F)
species 1 0.6488906 0.0757032 2.115638 0.002
type 1 0.5615418 0.0655126 1.830847 0.008
Residual 24 7.3610760 0.8587842 NA NA
Total 26 8.5715084 1.0000000 NA NA
                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.3294975 0.9684059 0.06064512   0.499      1.000    
2   Control vs Hot_control  1 0.7123847 2.3071902 0.11949901   0.003      0.009   *
3 Treatment vs Hot_control  1 0.4044284 1.3387544 0.07721168   0.054      0.162    
5.5.4.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.006755 0.0033773 0.4095    999   0.69
Residuals 24 0.197956 0.0082482                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.30200     0.957
Hot_control 0.30023                 0.565
Treatment   0.95485     0.53193          
Df SumOfSqs R2 F Pr(>F)
species 1 0.7984743 0.104299 3.065171 0.001
type 1 0.6051778 0.079050 2.323148 0.006
Residual 24 6.2519772 0.816651 NA NA
Total 26 7.6556293 1.000000 NA NA
                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.2860314 0.9858296 0.06166897   0.427      1.000    
2   Control vs Hot_control  1 0.8705308 3.2919398 0.16222894   0.003      0.009   *
3 Treatment vs Hot_control  1 0.4377241 1.6307998 0.09249721   0.037      0.111    
5.5.4.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00404 0.0020184 0.1345    999  0.909
Residuals 24 0.36024 0.0150100                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.84200     0.693
Hot_control 0.81332                 0.812
Treatment   0.63160     0.76897          
Df SumOfSqs R2 F Pr(>F)
species 1 0.0700374 0.0635900 1.6594454 0.15
type 1 0.0184254 0.0167292 0.4365646 0.80
Residual 24 1.0129278 0.9196808 NA NA
Total 26 1.1013906 1.0000000 NA NA
                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.01869516 0.4206549 0.02727867   0.781      1.000    
2   Control vs Hot_control  1 0.07762661 2.2806818 0.11828844   0.036      0.108    
3 Treatment vs Hot_control  1 0.03350315 0.6872008 0.04118131   0.668      1.000    
5.5.4.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01287 0.0064374 0.4319    999   0.66
Residuals 24 0.35771 0.0149044                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.38600     0.758
Hot_control 0.38202                 0.598
Treatment   0.75152     0.59510          
Df SumOfSqs R2 F Pr(>F)
species 1 -0.0048931 -0.0058792 -0.1628739 0.826
type 1 0.1161466 0.1395538 3.8660898 0.100
Residual 24 0.7210172 0.8663254 NA NA
Total 26 0.8322707 1.0000000 NA NA
                     pairs Df  SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.14708071 6.531088 0.30333294   0.043      0.129    
2   Control vs Hot_control  1 0.03067531 1.027084 0.05697451   0.341      1.000    
3 Treatment vs Hot_control  1 0.04199233 1.256701 0.07282392   0.301      0.903    
beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

5.5.4.3 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

5.5.4.4 Number of samples used

[1] 28
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)
5.5.4.4.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.001996 0.0009978 0.2043    999  0.825
Residuals 25 0.122093 0.0048837                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.70400     0.805
Hot_control 0.67713                 0.609
Treatment   0.79246     0.58926          
Df SumOfSqs R2 F Pr(>F)
type 2 1.542992 0.1948479 3.025016 0.001
Residual 25 6.375966 0.8051521 NA NA
Total 27 7.918958 1.0000000 NA NA
                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.2206497 0.7906964 0.04709134   0.862      1.000    
2   Control vs Hot_control  1 0.7211095 2.7164967 0.13777786   0.001      0.003   *
3 Treatment vs Hot_control  1 0.7163430 2.6326288 0.13409456   0.001      0.003   *
5.5.4.4.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.007233 0.0036166 0.6989    999  0.497
Residuals 25 0.129365 0.0051746                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.51500     0.645
Hot_control 0.49241                 0.283
Treatment   0.65999     0.28046          
Df SumOfSqs R2 F Pr(>F)
type 2 1.962645 0.2556667 4.293552 0.001
Residual 25 5.713932 0.7443333 NA NA
Total 27 7.676577 1.0000000 NA NA
                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.1081320 0.3901033 0.02380115   0.998      1.000    
2   Control vs Hot_control  1 0.7336222 2.8624966 0.14411565   0.001      0.003   *
3 Treatment vs Hot_control  1 0.7084005 2.6970385 0.13692609   0.001      0.003   *
5.5.4.4.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000406 0.0002032 0.0494    999  0.963
Residuals 25 0.102865 0.0041146                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.94000     0.831
Hot_control 0.93747                 0.787
Treatment   0.84169     0.75551          
Df SumOfSqs R2 F Pr(>F)
type 2 0.1793566 0.2173098 3.470559 0.001
Residual 25 0.6459931 0.7826902 NA NA
Total 27 0.8253496 1.0000000 NA NA
                     pairs Df   SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.005034708 0.178068 0.01100675   0.979      1.000    
2   Control vs Hot_control  1 0.099703372 3.519604 0.17152397   0.003      0.009   *
3 Treatment vs Hot_control  1 0.079909458 2.912000 0.14624347   0.004      0.012   .
5.5.4.4.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00836 0.0041823 0.2225    999  0.844
Residuals 25 0.46991 0.0187964                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.58600     0.679
Hot_control 0.53209                 0.855
Treatment   0.60231     0.82567          
Df SumOfSqs R2 F Pr(>F)
type 2 0.0031465 0.0045388 0.0569937 0.88
Residual 25 0.6900904 0.9954612 NA NA
Total 27 0.6932368 1.0000000 NA NA
                     pairs Df   SumsOfSqs   F.Model          R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.010979753 0.6382789 0.038362075   0.424          1    
2   Control vs Hot_control  1 0.003756213 0.1351199 0.007885556   0.702          1    
3 Treatment vs Hot_control  1 0.019226092 0.5508465 0.031385752   0.497          1    
beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")